%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% Assembly
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% define the ordering in the VAR
inflpos = 1;
y1pos = 2;
ysprpos = 3;
gdppos = 4;
divgrpos = 5;
pdpos = 6;
dtpos = 7;
coint_tax = 8;
dgpos = 9;
coint_spending = 10;
ddpos = 11;
coint_debt = 12;

%% raw data
clear Xraw;
Xraw(:, inflpos) = inflation;
Xraw(:, gdppos) = rpcgdpgr;
Xraw(:, y1pos) = ynom1;
Xraw(:, ysprpos) = yieldspr;
Xraw(:, pdpos) = pdm;
Xraw(:, divgrpos) = divgrm;
Xraw(:, dtpos) = deltalogtau;
Xraw(:, coint_tax) = log(taxrevgdp);
Xraw(:, dgpos) = deltalogg;
Xraw(:, coint_spending) = log(spendgdp);
Xraw(:, ddpos) = deltalogd;
Xraw(:, coint_debt) = logdts;

% VAR elements
infl = inflation - pi0;
x = rpcgdpgr - x0;
yield1 = ynom1 - y0nom_1;
nomyspr = yieldspr - yspr0;
dg = deltalogg - g0;
dt = deltalogtau - tau0;
pd = pdm - A0m;
% divgr = divgrm - mu_m;

%% trend adjustment

X2(:, inflpos) = infl;
X2(:, y1pos) = yield1;
X2(:, ysprpos) = nomyspr;
X2(:, gdppos) = x;
X2(:, divgrpos) = divgrm - mean(divgrm);
X2(:, pdpos) = pd;
X2(:, dtpos) = dt;
X2(:, coint_tax) = log(tts) - mean(log(tts));
X2(:, dgpos) = dg;
X2(:, coint_spending) = log(gts) - mean(log(gts));
X2(:, ddpos) = deltalogd_demean;
X2(:, coint_debt) = logdts;

N = cols(X2);
I = eye(N);
I_pi = I(:, inflpos);
I_gdp = I(:, gdppos);
I_y1 = I(:, y1pos);
I_yspr = I(:, ysprpos);
I_pdm = I(:, pdpos);
I_divgrm = I(:, divgrpos);
I_dt = I(:, dtpos);
I_dg = I(:, dgpos);
I_coint_tax = I(:, coint_tax);
I_coint_spending = I(:, coint_spending);
I_coint_debt = I(:, coint_debt);

z_var = X2(2:end, :);
z_varlag = X2(1:end - 1, :);
Y = z_var;
X = z_varlag;

% actual, non-detrended series of tax and spending
X2t = X2;
X2t(:, coint_tax) = log(taxrevgdp) - mean(log(taxrevgdp));
X2t(:, coint_spending) = log(spendgdp) - mean(log(spendgdp));

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% Estimate Psi and Sig
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

Psi = zeros(N, N);
Psi_se = zeros(N, N);

for i = [inflpos y1pos ysprpos gdppos divgrpos pdpos dtpos dgpos ddpos]
    regr = ols(Y(:, i), [ones(T - 1, 1), X(:, :)]);
    c(i) = regr.beta(1);
    c_se(i) = regr.bstd(1);
    Psi(i, :) = regr.beta(2:end)';
    Psi_se(i, :) = regr.bstd(2:end)';
    R2(i) = regr.rsqr;
    R2bar(i) = regr.rbar;
    eps(:, i) = Y(:, i) - c(i) - X * Psi(i, :)';

    clear regr;
end

Psi(coint_tax, :) = Psi(dtpos, :);
Psi(coint_spending, :) = Psi(dgpos, :);
Psi(coint_debt, :) = Psi(ddpos, :);
Psi(coint_tax, coint_tax) = Psi(coint_tax, coint_tax) + 1;
Psi(coint_spending, coint_spending) = Psi(coint_spending, coint_spending) + 1;

Psi(coint_debt, coint_debt) = Psi(coint_debt, coint_debt) + 1;

Psi_se(coint_tax, :) = Psi_se(dtpos, :);
Psi_se(coint_spending, :) = Psi_se(dgpos, :);

Psi_se(coint_debt, :) = Psi_se(ddpos, :);

Sigma = cov(eps(:, [inflpos y1pos ysprpos gdppos divgrpos pdpos dtpos dgpos]));
tmp = chol(Sigma, 'lower');
Sig = zeros(N, N);
Sig([inflpos y1pos ysprpos gdppos divgrpos], [inflpos y1pos ysprpos gdppos divgrpos]) = ...
    tmp(1:5, 1:5);
Sig(pdpos, [inflpos y1pos ysprpos gdppos divgrpos pdpos]) = tmp(6, 1:6);
Sig(dtpos, [inflpos y1pos ysprpos gdppos divgrpos pdpos dtpos]) = tmp(7, 1:7);
Sig(dgpos, [inflpos y1pos ysprpos gdppos divgrpos pdpos dtpos dgpos]) = tmp(8, 1:8);

Sig(coint_tax, :) = Sig(dtpos, :);
Sig(coint_spending, :) = Sig(dgpos, :);

eps = [eps(:, [inflpos y1pos ysprpos gdppos divgrpos]), ...
                eps(:, pdpos), eps(:, dtpos), eps(:, dtpos), ...
                eps(:, dgpos), eps(:, dgpos), eps(:, ddpos), eps(:, ddpos)];
Sigma = Sig * Sig';

tstat = Psi ./ Psi_se;

max(abs(eig(Psi)))

eps2 = eps;

tau0 = 0;
g0 = 0;
